#include <stdio.h>
#include <stdlib.h>
#include <limits.h>
#include <pthread.h> // Colocar -pthread no gcc
#include <time.h>
#include <gmp.h> // Tem que colocar -lgmp na hora de compilar

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

#define N_THREADS 2*N

void init_mpf_b(mpf_t rop);
void media_simples(mpf_t rop, mpf_t a, mpf_t b);
void media_geometrica(mpf_t rop, mpf_t a, mpf_t b); 
void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a);
void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t);

pthread_mutex_t MUTEX;
pthread_cond_t isItOk;

mpf_t a;
mpf_t b;
mpf_t t;
mpf_t x;
mpf_t y;
mpf_t pi;

void *ThreadA(void *tid){
    int rc;

    pthread_mutex_lock(&MUTEX);
            
    media_simples(a, a, b);     //a = (a+b)/2.0;

    pthread_mutex_unlock(&MUTEX);
    pthread_cond_signal(&isItOk);  
   
    pthread_cond_wait(&isItOk, &MUTEX);

    media_geometrica(b, b, y);  //b = sqrt(b*y); 
    mpf_set(y, a);              //y = a;

    pthread_mutex_unlock(&MUTEX);
    
    pthread_exit(NULL);
}

void *ThreadB(void *tid){


    equacao_t(t, x, y, a);      //t = t -   x * pow((y-a),2); 
    mpf_mul_ui(x, x, 2);        //x *= 2; 
   
    pthread_cond_signal(&isItOk);  
    pthread_mutex_unlock(&MUTEX);
    
pthread_exit(NULL);
}

int main(int argc, char* argv[]){

    unsigned long i, j=0, p=0, q=0;

    pthread_t threads[N_THREADS];
//    char alive[N_THREADS];
    pthread_attr_t attr;
    void *exit_status;

    time_t t_zero, t_final;
    double delta_t;

    // Marca o tempo
    time(&t_zero);

//    memset(alive, '\0', sizeof(char)*N_THREADS);

    // Colocar Pthreads como Joinable    
    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
    // inicializar Mutexes
    pthread_cond_init(&isItOk, NULL);
    pthread_mutex_init(&MUTEX, NULL);

    // Float com precisao arbitraria
    mpf_set_default_prec(MANTISSA);

    mpf_init(a);
    mpf_init(b);
    mpf_init(t);
    mpf_init(x);
    mpf_init(y);
    mpf_init(pi);

    //  a = 1.0;
    //  b = (1.0/(sqrt(2)));
    //  t = (1.0/4.0);
    //  x = 1.0;
    //  pi=3.14, y=1.0; // Independnetes

    mpf_set_d(t, (1.0/4.0));
    mpf_set_ui(a, 1);
    mpf_set_ui(y, 1);
    mpf_set_ui(x, 1);
    init_mpf_b(b);

    /*
     *  A(n+1) = A(n)+B(n)/2
     *  B(n+1) = sqrt(A(n)+B(n))
     *  T(n+1) = T(n) - P(n)(A(n) - A(n+1))^2
     *  P(n+1) = 2*P(n)
     *  pi = ((A(n)+B(n))^2)/4T(n)
     * 
     * */

    for(i=0; i<N; i++){
       
        p=j;
        pthread_create(&threads[j], &attr, ThreadA, (void*) q);    
        j++; q++;

        pthread_cond_wait(&isItOk, &MUTEX);
        pthread_create(&threads[j], &attr, ThreadB, (void*) q);
        j++; q++;
        
        for(j=p;j<q; j++){
            pthread_join(threads[j], NULL);
        }
    }
    

    equacao_pi(pi, a, b, t);    //pi = pow((a+b),2)/(4.0*t);
    
    // Limpa a Caquinha
    mpf_clear(a);
    mpf_clear(b);
    mpf_clear(y);
    mpf_clear(t);

    mpf_out_str(NULL, BASE, 0, pi);
    mpf_clear(pi);
    
    pthread_attr_destroy(&attr);
    pthread_mutex_destroy(&MUTEX);
    pthread_cond_destroy(&isItOk);

    time(&t_final);
    delta_t = difftime(t_final, t_zero);
    printf("\nTempo = %.2lf segundos\n", delta_t);

    pthread_exit(NULL);   
    
    return 0;
}

void init_mpf_b(mpf_t rop){
    
    mpf_set_default_prec(MANTISSA);
    mpf_t x, y;

    mpf_init(x);
    mpf_init(y);

    // x=1
    mpf_set_ui(x, 1);
    // y=sqtr(2)
    mpf_sqrt_ui(y, 2);
    // x=1/sqrt(2)
    mpf_div(x, x, y);

    // return
    mpf_set(rop, x);

    mpf_clear(x);
    mpf_clear(y);
}

void media_simples(mpf_t rop, mpf_t a, mpf_t b){
    mpf_t aux;  
    mpf_init(aux);
    //aux=a+b;
    mpf_add(aux, a, b);
    //rop = aux/2
    mpf_div_ui(rop, aux, 2);
    mpf_clear(aux);
}

void media_geometrica(mpf_t rop, mpf_t a, mpf_t b){ 
    mpf_t aux;
    mpf_init(aux);
    //aux=a+b;
    mpf_mul(aux, a, b);
    //rop = sqrt(aux)
    mpf_sqrt(rop, aux);
    mpf_clear(aux);
}

void equacao_t(mpf_t t, mpf_t x, mpf_t y, mpf_t a){
    mpf_t aux;
    mpf_init(aux);
    // t = t - x * pow((y-a),2);
    mpf_sub(aux, y, a);
    mpf_pow_ui(aux, aux, 2);
    mpf_mul(aux, aux, x);
    mpf_sub(t, t, aux);
    mpf_clear(aux);
}

void equacao_pi(mpf_t pi, mpf_t a, mpf_t b, mpf_t t){
    mpf_t aux;
    mpf_init(aux);
    //pi = pow((a+b),2)/(4.0*t);
    mpf_add(aux, a, b);
    mpf_pow_ui(aux, aux, 2);
    mpf_div(aux, aux, t);
    mpf_div_ui(pi, aux, 4);
    mpf_clear(aux);
}
